Add gene names to the result using the same gtf file

res.Salmon.basic$gene <- gtf.data$gene_name[match(rownames(res.Salmon.basic), gtf.data$gene_id)]
res.Salmon.basic$ensembl <- gsub('\\..+$', '', rownames(res.Salmon.basic))
res.Salmon.basic$entrez <- AnnotationDbi::mapIds(org.Hs.eg.db::org.Hs.eg.db, keys=res.Salmon.basic$ensembl, column="ENTREZID", keytype="ENSEMBL", multiVals = "first")

## 'select()' returned 1:many mapping between keys and columns

Log fold change (LFC) shrinkage

Shrinkage of effect size (LFC estimates) is useful for visualization and ranking of genes. ‘apeglm’ and ‘ashr’ are preferred method over ‘normal’. Please read more about each method from the original papers (cited below) as well as from the Alternative shrinkage estimators and Extended section on shrinkage estimators. As for deciding between apeglm vs. ashr, see this thread and this thread.

## LFC shrinkage methods : type = c("apeglm", "ashr", "normal")
resLFC.Salmon.basic.apeglm <- DESeq2::lfcShrink(dds.Salmon.basic.filtered, res = res.Salmon.basic, type = "apeglm", 
                                                coef = "condition_CC_CR_vs_CC")

## Warning in nbinomGLM(x = x, Y = YNZ, size = size, weights = weightsNZ, offset =
## offsetNZ, : the line search routine failed, unable to sufficiently decrease the
## function value

resLFC.Salmon.basic.ashr <- DESeq2::lfcShrink(dds.Salmon.basic.filtered, res = res.Salmon.basic, type = "ashr")
resLFC.Salmon.basic.normal <- DESeq2::lfcShrink(dds.Salmon.basic.filtered, res = res.Salmon.basic, type = "normal",
                                                coef = "condition_CC_CR_vs_CC")

## reorder by adjusted p-value
resLFC.Salmon.basic.apeglm <- resLFC.Salmon.basic.apeglm[order(resLFC.Salmon.basic.apeglm$padj),]
resLFC.Salmon.basic.ashr <- resLFC.Salmon.basic.ashr[order(resLFC.Salmon.basic.ashr$padj),]
resLFC.Salmon.basic.normal <- resLFC.Salmon.basic.normal[order(resLFC.Salmon.basic.normal$padj),]

## MA plot before/after shrinkage
Plot.MA <- function(ma.data, baseMeanCutoff = 10, FDR = 0.05, title = "", colSig = "black") {

  ma <- data.frame(ma.data)
  ma <- ma %>% dplyr::mutate(isDE = case_when(
    abs(log2FoldChange) > 0.5 & baseMean >= baseMeanCutoff & padj < FDR ~ TRUE,
    TRUE ~ FALSE),
    isDE = as.factor(isDE))
  
  ma <- ma %>% dplyr::arrange(isDE)
  
  ma.plot <- ma %>% data.frame() %>% ggplot2::ggplot(aes(x=baseMean, y=log2FoldChange, color=isDE, shape=isDE, fill=isDE)) +
    geom_point(aes(shape = isDE, color=isDE), size = 3) +
    scale_x_continuous(trans = "log10") +
    scale_color_manual(values = c("grey80", colSig)) +
    scale_shape_manual(values= c(19,16)) +
    geom_vline(xintercept = 50, color = 'grey50', linetype = "dashed") +
    geom_hline(yintercept = c(-.5, .5), color = 'grey50', linetype = "dashed") +
    Theme.YL + ggtitle(label = title)
    
  plot(ma.plot)
  
}

Plot.MA(res.Salmon.basic, title = paste("MA plot (before shrinkage)", 
              "min:", round(range(res.Salmon.basic$log2FoldChange),2)[1], 
              "max:", round(range(res.Salmon.basic$log2FoldChange),2)[2]), 
        colSig = color$Extra[1])

Plot.MA(resLFC.Salmon.basic.apeglm, title = paste("MA plot (after apeglm shrinkage)", 
                                           "min:", round(range(resLFC.Salmon.basic.apeglm$log2FoldChange),2)[1], 
                                           "max:", round(range(resLFC.Salmon.basic.apeglm$log2FoldChange),2)[2]), 
        colSig = color$Extra[2])

Plot.MA(resLFC.Salmon.basic.ashr, title = paste("MA plot (after ashr shrinkage)", 
                                                "min:", round(range(resLFC.Salmon.basic.ashr$log2FoldChange),2)[1], 
                                                "max:", round(range(resLFC.Salmon.basic.ashr$log2FoldChange),2)[2]), 
        colSig = color$Extra[3])

Plot.MA(resLFC.Salmon.basic.normal, title = paste("MA plot (after normal shrinkage)", 
                                                  "min:", round(range(resLFC.Salmon.basic.normal$log2FoldChange),2)[1], 
                                                  "max:", round(range(resLFC.Salmon.basic.normal$log2FoldChange),2)[2]), 
        colSig = color$Extra[4])

Save DEA result

## choose one shrinkage method
resLFC.Salmon.basic <- resLFC.Salmon.basic.normal
rm(resLFC.Salmon.basic.apeglm, resLFC.Salmon.basic.ashr)

write.table(resLFC.Salmon.basic, file = paste0(outdir,"DEA_all_genes_shrunkLFC.tsv"), sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)

ipa <- resLFC.Salmon.basic[,c("gene","ensembl","entrez", "baseMean", "log2FoldChange","pvalue", "padj")]
write.table(ipa, file = paste0(outdir,"DEA_all_genes_IPA.txt"), sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)

Volcano Plot

A volcano plot allows quick visual identification of genes with large fold changes that are also statistically significant. These may be the most biologically significant genes.

mean.cutoff = 10
lfc.cutoff = 0.5 
threshold <- resLFC.Salmon.basic$padj < FDR & 
  resLFC.Salmon.basic$baseMean > mean.cutoff & 
  abs(resLFC.Salmon.basic$log2FoldChange) > lfc.cutoff

df <- data.frame(resLFC.Salmon.basic) 
df$Significant <- threshold
df$genelabels <- rownames(df) %in% rownames(df)[df$Significant == TRUE][1:10]

p <- ggplot(df) +
  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = Significant, alpha = Significant, size = Significant)) +
  scale_color_manual(values=color$Extra) +
  geom_label_repel(aes(x = log2FoldChange, y = -log10(padj), label = ifelse(genelabels == T, df$gene,""))) +
  labs(title = "Volcano Plot (CC_CR vs. CC)", 
       subtitle = paste0("10 significant genes are labelled. LFC cutoff = ", lfc.cutoff, 
                         ", p-value cutoff = ", FDR, ", baseMean cutoff = ", mean.cutoff)) + 
  xlab("log2 fold change") +  ylab("-log10 adjusted p-value") + 
  geom_vline(xintercept = c(lfc.cutoff, -lfc.cutoff),  colour = "grey50", linetype = 'dashed') + 
  geom_hline(yintercept = FDR, colour = "grey50", linetype = 'dashed') + 
  Theme.YL + scale_size_manual(values = c(1,2))
plot(p)

## Warning: Using alpha for a discrete variable is not advised.

## Warning: Use of `df$gene` is discouraged.
## ℹ Use `gene` instead.

## png 
##   2

PlotCounts

Check out counts for a gene of your interest. Normalized counts plus a pseudocount of 0.5 are shown by default.

gene.interest.name =  "MIR100HG" # high variability within CC_CR group 
gene.interest = rownames(resLFC.Salmon.basic[resLFC.Salmon.basic$gene == gene.interest.name,]) 

d <- plotCounts(dds.Salmon.basic.filtered, gene = gene.interest, intgroup = "condition", returnData = TRUE)
d$name <- rownames(d)
d$tech_rep <- SampleSheet$tech_rep

ggplot(d, aes(x=condition, y=count, color=tech_rep)) + scale_color_manual(values=color$Replicate) + scale_fill_manual(values=color$Condition) +
  geom_point(position=position_jitter(w=0.1,h=0)) + geom_boxplot(aes(group = condition, fill= condition),alpha = 0.3) + #geom_text_repel(aes(label = name)) +
  ggtitle(paste(gene.interest.name, "")) + Theme.YL

gene.interest.name =  "RPL5" # low variability within groups
gene.interest = rownames(resLFC.Salmon.basic[resLFC.Salmon.basic$gene == gene.interest.name,]) 

d <- plotCounts(dds.Salmon.basic.filtered, gene = gene.interest, intgroup = "condition", returnData = TRUE)
d$name <- rownames(d)
d$tech_rep <- SampleSheet$tech_rep

ggplot(d, aes(x=condition, y=count, color=tech_rep)) + scale_color_manual(values=color$Replicate) + scale_fill_manual(values=color$Condition) +
  geom_point(position=position_jitter(w=0.1,h=0)) + geom_boxplot(aes(group = condition, fill= condition),alpha = 0.3) + #geom_text_repel(aes(label = name)) +
  ggtitle(paste(gene.interest.name, "")) + Theme.YL

Data Transformation

For visualization or clustering, it might be useful to work with transformed versions of the count data. The regularized log (rlog) function transforms the count data to the log2 scale in a way which minimizes differences between samples for rows with small counts, and which normalizes with respect to library size. The rlog transformation produces a similar variance stabilizing effect as Variance Stabilizing Transformation (vst), though rlog is more robust in the case when the size factors vary widely.

TransformCounts <- function(ddsMatrix, TransformationType = "rlog", blind = TRUE, fitType = "parametric") {
  ## user can choose between rlog and vst as their data transformation method
  ## check if the ddsMatrix contains already calculated dispersion values or not
    if (sum(is.na(SummarizedExperiment::rowData(ddsMatrix)[,"dispersion"])) == 0) {
      transformedCounts <- switch(
        TransformationType,
        "rlog" = DESeq2::rlog(ddsMatrix, fitType = fitType, blind = blind),
        "vst" = DESeq2::varianceStabilizingTransformation(ddsMatrix, fitType = fitType, blind = blind)) 
      } else {
        transformedCounts <- switch(
          TransformationType,
          "rlog" = DESeq2::rlog(estimateDispersions(estimateSizeFactors(ddsMatrix)), fitType = fitType, blind = blind), 
          "vst" = DESeq2::varianceStabilizingTransformation(estimateDispersions(estimateSizeFactors(ddsMatrix)), fitType = fitType, blind = blind))
          }
  hist(assay(ddsMatrix))
  hist(assay(transformedCounts))
  return(transformedCounts)
}

## all genes : N = 34378 
rlog.dds.Salmon.basic <- TransformCounts(dds.Salmon.basic.filtered, "rlog")

## Criteria : baseMean > 10 AND adjusted p-value < 0.05 AND |log2FoldChange| > 0.5
# replace NA values to FC:0 and padj:1, then select significant genes
res.Salmon.basic.Filter <- res.Salmon.basic %>% data.frame() %>%
  dplyr::mutate_at(vars(contains("log2FoldChange")), ~replace(., is.na(.), 0)) %>%
  dplyr::mutate_at(vars(contains("padj")), ~replace(., is.na(.), 1)) %>%
  dplyr::filter(baseMean > mean.cutoff, padj < FDR, abs(log2FoldChange) > lfc.cutoff)
res.Salmon.basic.Filter <- res.Salmon.basic.Filter[order(res.Salmon.basic.Filter$padj),]

## only the significant genes : N = 831 
rlog.dds.Salmon.basic.sig <- TransformCounts(dds.Salmon.basic.filtered[rownames(res.Salmon.basic.Filter),], "rlog") 

dds.Salmon.basic.filtered[rownames(res.Salmon.basic.Filter),]

## class: DESeqDataSet 
## dim: 1008 24 
## metadata(1): version
## assays(8): counts avgTxLength ... replaceCounts replaceCooks
## rownames(1008): ENSG00000122406.14 ENSG00000117523.18 ...
##   ENSG00000186567.14 ENSG00000289463.1
## rowData names(23): baseMean baseVar ... maxCooks replace
## colnames(24): CC_rep1_s1 CC_rep1_s2 ... CC_CR_rep3_s3 CC_CR_rep3_s4
## colData names(4): condition bio_rep tech_rep replaceable

Visualizing clusters with heatmap and dendogram

Using the count matrix with all genes, rlog transformation had been performed, and Euclidean distances among the samples are visualized as a heatmap and hierarchical clustering is performed using complete linkage method. Learn more about how to use pheatmap package: click here

## calculate distance 
sampleDists <- dist(t(assay(rlog.dds.Salmon.basic)), method = "euclidean")
sampleDistMatrix <- as.matrix(sampleDists)

## annotation color
mat_col <- data.frame(Condition = rlog.dds.Salmon.basic$condition, Replicate = rlog.dds.Salmon.basic$tech_rep)
rownames(mat_col) <- colnames(sampleDistMatrix)
mat_colors <- list(Condition = color$Condition, Replicate = color$Replicate)
names(mat_colors$Condition) <- unique(mat_col$Condition)
names(mat_colors$Replicate) <- unique(mat_col$Replicate)

## heatmap
p <- ggplotify::as.ggplot(pheatmap(mat = sampleDistMatrix,
         clustering_distance_rows=sampleDists, clustering_distance_cols=sampleDists,
         clustering_method = "complete",
         color = colorRampPalette(rev(brewer.pal(n = 9, name =  "YlGnBu")))(255),
         border_color = NA, cutree_cols=2, show_rownames = FALSE, fontsize = 15, #fontfamily = "Sans",
         width = ncol(sampleDistMatrix)*unit(7, "mm"), height = nrow(sampleDistMatrix)*unit(7, "mm"),
         cellwidth = 20, cellheight = 16,
         annotation_col = mat_col, annotation_colors = mat_colors, 
         main = "Clustering heatmap using all genes"))

## png 
##   2

Same as the plot above, but using the count matrix with only the significant genes.

## calculate distance 
sampleDists <- dist(t(assay(rlog.dds.Salmon.basic.sig)), method = "euclidean")
sampleDistMatrix <- as.matrix(sampleDists)

## annotation color
mat_col <- data.frame(Condition = rlog.dds.Salmon.basic.sig$condition, Replicate = rlog.dds.Salmon.basic.sig$tech_rep)
rownames(mat_col) <- colnames(sampleDistMatrix)
mat_colors <- list(Condition = color$Condition, Replicate = color$Replicate)
names(mat_colors$Condition) <- unique(mat_col$Condition)
names(mat_colors$Replicate) <- unique(mat_col$Replicate)


## heatmap
p <- ggplotify::as.ggplot(pheatmap(mat = sampleDistMatrix,
         clustering_distance_rows=sampleDists, clustering_distance_cols=sampleDists,
         clustering_method = "complete",
         color = colorRampPalette(rev(brewer.pal(n = 9, name =  "YlGnBu")))(255),
         border_color = NA, cutree_cols=2, show_rownames = FALSE, fontsize = 15, #fontfamily = "Sans",
         width = ncol(sampleDistMatrix)*unit(7, "mm"), height = nrow(sampleDistMatrix)*unit(7, "mm"),
         cellwidth = 20, cellheight = 16,
         annotation_col = mat_col, annotation_colors = mat_colors, 
         main = "Clustering heatmap using Significant genes"))

## png 
##   2

PCA and t-SNE plot

To make the large dataset of 34395 rows and 6 columns, dimensionality reduction techniques like PCA (Principal component analysis) and t-SNE (t-distributed stochastic neighbor embedding) are utilized for data visualization. We can utilize the clustering results as quality control measure of the experiment.

Scaling (z-transformation)

z_Transform <- function(rlog_dds) {
  normalized.counts = SummarizedExperiment::assay(rlog_dds)
  ## gene-wise z-transformation
  m = apply(normalized.counts, 1, mean)
  s = apply(normalized.counts, 1, sd)
  z = data.frame((normalized.counts - m)/s)
  ## add gene names
  # gtf.data <- GenomicTools.fileHandler::importGTF(gtf.path, level="gene", features=c("gene_id", "gene_name", "gene_type"))
  z$gene = gtf.data$gene_name[match(rownames(rlog_dds), gtf.data$gene_id)]
  z$ensembl <- gsub('\\..+$', '', rownames(rlog_dds))
  z$entrez <- AnnotationDbi::mapIds(org.Hs.eg.db::org.Hs.eg.db, keys=z$ensembl, column="ENTREZID", keytype="ENSEMBL", multiVals = "first")
  return(z)
}


z.counts <- z_Transform(rlog.dds.Salmon.basic)

## 'select()' returned 1:many mapping between keys and columns

write.csv(z.counts, file = paste0(outdir,"gene-wise-z-transformation.csv"))

Recreation of Figure 2a

Heatmap of the top 50 differentially expressed transcripts in CC-CR versus CC from three independent 3D culture experiments. Gene expression values are gene-wise z-transformed and are colored red for high abundance and blue for low abundance, as indicated in the scale bar.

## Top 50 genes from the original paper (fig 2a)
top50 = c("MIR100HG","HS3ST5","A2M","ZNF804A", "ARHGAP29", "SESN3", "SATB1", "SERPINE2", "AKAP12", "EPHB6", "GAS6", "IGFBP6", "CUEDC1", "CGN", "MESDC1", "VIL1", "SNTB1", "RP11-395P17.3", "ALPK3", "LOXL1", "CAPS", "FYN", "H19", "RASGRF1", "CDK6", "PKP1", "DYSF", "HLA-DRA", "ENC1", "QPCT", "IGF2BP2", "AC007952.5", "SPDYC", "HOXD10", "CIB2", "CFTR", "HOXD11", "TEX30", "LINC01018", "SERTAD4", "DKK3", 'GATA6', 'GSTM4', "CPVL", 'DKK1', 'EPB41L3', 'ZNK608', "CRABP1", "SERP2", "OLFM1") 

z.counts.top50 <- z.counts[z.counts$gene %in% top50,]
z.counts.top50 <- z.counts.top50 %>% arrange(factor(gene, levels = top50))
rownames(z.counts.top50) <- z.counts.top50$gene
z.counts.top50$gene <- NULL
z.counts.top50$ensembl <- NULL
z.counts.top50$entrez <- NULL

## Heatmap using gene-wise z-transformed gene expression values
mat_col <- data.frame(sample = rownames(coldata), 
                      Condition = coldata$condition)
rownames(mat_col) <- mat_col$sample
mat_col$sample <- NULL

p <- pheatmap(z.counts.top50, cluster_rows = FALSE,
         border_color = NA, fontsize = 15, 
         width = ncol(z.counts.top50)*unit(7, "mm"), height = nrow(z.counts.top50)*unit(7, "mm"),
         cellwidth = 20, cellheight = 16,
         annotation_col = mat_col, annotation_colors = mat_colors,
         main = "Top 50 differentially expressed genes in CC_CR vs. CC")

PlotCounts of MIR100HG gene

gene.interest.name = "MIR100HG"
gene.interest = rownames(resLFC.Salmon.basic[resLFC.Salmon.basic$gene == gene.interest.name,]) 

d <- plotCounts(dds.Salmon.basic.filtered, gene = gene.interest, intgroup = "condition", returnData = TRUE)
d$name <- rownames(d)
d$tech_rep <- SampleSheet$tech_rep

ggplot(d, aes(x=condition, y=count, color=tech_rep)) + scale_color_manual(values=color$Replicate) + scale_fill_manual(values=color$Condition) +
  geom_point(position=position_jitter(w=0.1,h=0)) + geom_boxplot(aes(group = condition, fill= condition),alpha = 0.3) + #geom_text_repel(aes(label = name)) +
  ggtitle(paste(gene.interest.name, "")) + Theme.YL

## png 
##   2

Ranking LFCs

## Remove NAs 
gseaDat <- resLFC.Salmon.basic[!is.na(resLFC.Salmon.basic$entrez),] 

## Remove duplicated gene names -- 1) significant? 2) highest baseMean
length(gseaDat$entrez[duplicated(gseaDat$entrez)])

## [1] 83

## 1) significant? 
gseaDat <- subset(gseaDat, abs(log2FoldChange) > lfc.cutoff & padj < FDR)

## No duplicates
length(gseaDat$entrez[duplicated(gseaDat$entrez)])

## [1] 0

## Rank all genes based on their shrunken fold change 
geneList <- gseaDat$log2FoldChange
names(geneList) <- gseaDat$entrez

## Plot the ranked fold changes in descending order
geneList <- sort(geneList, decreasing = T)
p <- barplot(geneList, main = "Ranked fold changes of all genes in descending order")

# pdf(paste0(outdir,"plot_rank_LFC.pdf"), width=10, height=8)
# plot(p)
# dev.off()

MSigDB & clusterProfiler

## Choose a category from MSigDB
MSigDB <- msigdbr::msigdbr(species = "Homo sapiens", category = "H")
MSigDB <- MSigDB %>% as.data.frame %>% dplyr::select(gs_name, entrez_gene)

set.seed(20250319)
gsea <- clusterProfiler::GSEA(geneList, TERM2GENE = MSigDB,
                                     minGSSize    = 10,
                                     maxGSSize    = 500,
                                     pvalueCutoff = 0.05,
                                     pAdjustMethod = "BH")

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

## preparing geneSet collections...

## GSEA analysis...

## leading edge analysis...

## done...

## significant gene set
gsea@result %>% dplyr::select(ID, NES, qvalue, setSize, rank) %>% dplyr::arrange(qvalue) %>% head()

##                                                                    ID      NES
## HALLMARK_INTERFERON_GAMMA_RESPONSE HALLMARK_INTERFERON_GAMMA_RESPONSE 2.683109
## HALLMARK_INTERFERON_ALPHA_RESPONSE HALLMARK_INTERFERON_ALPHA_RESPONSE 2.368676
##                                          qvalue setSize rank
## HALLMARK_INTERFERON_GAMMA_RESPONSE 0.0004538009      25  147
## HALLMARK_INTERFERON_ALPHA_RESPONSE 0.0017829376      19  335

openxlsx::write.xlsx(gsea@result %>% dplyr::arrange(qvalue), 
                     file = paste0(outdir, "/gsea.hallmark.xlsx"))

GeneRatio = count / setSize. ‘count’ is the number of genes that belong to a given gene-set, while ‘setSize’ is the total number of genes in the gene-set.

enrichplot::dotplot(gsea, showCategory = 10) + 
  ggplot2::scale_fill_gradient(low = "#132B43",  high = "#56B1F7", limits = c(0, 1), name = "qvalue") 

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

enrichplot::gseaplot(gsea, geneSetID = 1, title = gsea$Description[1])

enrichplot::gseaplot(gsea, geneSetID = 1, title = gsea$Description[2])

reactomePathways

## Pathway analysis
pathways <- fgsea::reactomePathways(names(geneList))  
fgseaRes <- fgsea::fgsea(pathways, geneList, maxSize = 500)

## Top 10 pathways
head(fgseaRes[order(padj, -abs(NES)), ], n=10)[,c(1,3,5)]

##                                          pathway      padj         ES
##                                           <char>     <num>      <num>
##  1:              Interferon alpha/beta signaling 0.3896538  0.5227749
##  2:  Antiviral mechanism by IFN-stimulated genes 0.3896538  0.5275529
##  3:                       Platelet degranulation 0.3896538  0.6107831
##  4: Response to elevated platelet cytosolic Ca2+ 0.3896538  0.6107831
##  5:                    ISG15 antiviral mechanism 0.3896538  0.5944044
##  6:                                 Gastrulation 0.3896538 -0.6505853
##  7:      Disorders of transmembrane transporters 0.3896538 -0.6820418
##  8:                  FOXO-mediated transcription 0.3896538  0.5555556
##  9:                         RHO GTPase Effectors 0.3896538 -0.6323468
## 10:             Crosslinking of collagen fibrils 0.3896538 -0.7974026

## Enrichment plot of the first pathway
pathway.interest = head(fgseaRes[order(padj, -abs(NES)), ], n=10)$pathway[1]
p <- fgsea::plotEnrichment(pathways[[pathway.interest]], geneList) + ggplot2::labs(title=pathway.interest) + Theme.YL

# pdf(paste0(outdir,"plot_pathway_", gsub('\\s|\\/', '_', pathway.interest), ".pdf"), width=10, height=8)
plot(p)

# dev.off()